Leçon 4: Python, calcul symbolique avec sympy.

Équations différentielles ordinaires.

La solution analytique des équations différentielles à condition initiale par exemple:

import sympy as sp
from IPython.display import *
%matplotlib inline
sp.init_printing(use_latex=True)
t,C1,C2=sp.symbols('t,C1,C2')
x=sp.Function('x')(t)
eq=sp.Eq(x.diff(t,2)-9.81)         # équation ballistique, trajectoire
display(eq)
xsol=sp.dsolve(eq,x).rhs
display(xsol)
constantes={'C1':0,'C2':0}      # les conditions initiales not mises ici dans le dictionnaire
                                # constantes qui a été défini. C1 est la position initiale 
                                # et C2 la vitesse initiale.
xsol=xsol.subs(C1,constantes['C1'])
xsol=xsol.subs(C2,constantes['C2'])
display(xsol)
display()

$$\frac{d^{2}}{d t^{2}} x{\left (t \right )} - 9.81 = 0$$

$$C_{1} + C_{2} t + 4.905 t^{2}$$

$$4.905 t^{2}$$

Voyons la solution d'une équation différentielle d'ordre 2 à conditions aux limites. Les 4 premières lignes sont utilisées comme dans la leçon précédente. (importer sympy et faire l'affichage de qualité en LaTeX)

import sympy as sp
from IPython.display import *
%matplotlib inline
sp.init_printing(use_latex=True)
x,N=sp.symbols('x,N')
T=sp.Function("T")(x)
eq=sp.Eq(T.diff(x,2)-N**2*T)                 # Transport Phenomena, 10.7, ailette droite
                                             # forme adimensionnelle
Sol=sp.dsolve(eq,T).rhs                      # Solution generale
display(Sol)

$$C_{1} e^{- N x} + C_{2} e^{N x}$$

On veut maintenant trouver les valeurs des constantes. On utilisera donc les conditions aux limites, T=1 en x=0 et dT/dx=0 en x=1.

                                             # Trouve C1 et C2 a partir des CL, 2 eq 2 inconnues
S=sp.solve([sp.Eq(Sol.subs(x,0),1),          # condition de Dirichlet 
            sp.Eq(Sol.diff(x).subs(x,1),0)]  # condition de Newman
            ,sp.symbols('C1,C2'))            # 
Final=Sol.subs(S)
display(Final)
final=sp.simplify(Final)
display(final)

$$\frac{e^{2 N} e^{- N x}}{e^{2 N} + 1} + \frac{e^{N x}}{e^{2 N} + 1}$$

$$\frac{e^{- N x}}{e^{2 N} + 1} \left(e^{2 N} + e^{2 N x}\right)$$

Autre exemple, un film dont la viscosité est variable s'écoulant sur un plan

import sympy as sp
from IPython.display import *
%matplotlib inline
sp.init_printing(use_latex=True)
rho,cosbeta,g,alpha,delta,mu0=sp.symbols('rho,cosbeta,g,alpha,delta,mu0',positive=True)
x=sp.symbols('x')
vz=sp.symbols('vz',function=True)
eq=mu0*sp.exp(-alpha*x/delta)*sp.Derivative(vz(x),x)+rho*g*cosbeta*x
display(eq)
vz=sp.dsolve(eq)
display(vz)

$$cosbeta g \rho x + \mu_{0} e^{- \frac{\alpha x}{\delta}} \frac{d}{d x} \operatorname{vz}{\left (x \right )}$$

$$\operatorname{vz}{\left (x \right )} = \frac{1}{\mu_{0}} \left(C_{1} - \frac{cosbeta x}{\alpha} \delta g \rho e^{\frac{\alpha x}{\delta}} + \frac{cosbeta g}{\alpha^{2}} \delta^{2} \rho e^{\frac{\alpha x}{\delta}}\right)$$

Je vous invite à explorer les nombreux documents disponibles sur l'utilisation de sympy, à commencer avec live.sympy.org , à voir absolument. La solution des équations différentielles de plusieurs types sont présentées.